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We study the quantum spin liquid phase in a variant of the Kitaev model where the bonds of the honeycomb 
lattice are distributed in a Kekule pattern. The system supports gapped and gapless Z 2 quantum spin liquids 
with interesting differences from the original Kitaev model, the most notable being a gapped Z 2 spin liquid on 
a Kagome lattice. Perturbing the exactly solvable model with antiferromagnetic Heisenberg perturbations, we 
find a magnetically ordered phase stabilized by a quantum ‘order by disorder’ mechanism, as well as an exotic 
continuous quantum phase transition between the topological spin liquid and this magnetically ordered phase. 
Using a combination of field theory and Monte-Carlo simulations, we find that the transition likely belongs to 
the 3D-XY x Z 2 universality class. 


I. INTRODUCTION 

Quantum spin liquids (QSL) represent prototypical con¬ 
densed matter phases whose description requires under¬ 
standing beyond the paradigm of spontaneous symmetry 
breaking.A QSL is a quantum paramagnet that can support 
quasiparticle excitations carrying quantum numbers which are 
fractions of the underlying microscopic degrees of freedom 
and hence are fundamentally different from random single 
spin flips of a thermal paramagnet or spin waves in a mag¬ 
netically ordered state.Systematic understanding of such 
phases, and phase transitions involving them, form an impor¬ 
tant area of current research in condensed matter physics. 

An important development in the understanding of QSLs 
came with the advent of exactly solvable spin Hamiltoni¬ 
ans where the ground state is a QSL and low energy exci¬ 
tations are indeed fractionalized. Lollowing the pioneering 
work of Kitaev,^’several such models are now known^^”^^ 
and their investigations have enhanced our understanding of 
QSL phases. These Kitaev models usually do not have spin 
rotation symmetry and the suggestion that some of them may 
be realized in 5d transition metal compounds (like Iridates), 
due to the presence of strong spin-orbit coupling, has lead to a 
plethora of interesting studies regarding their properties. 

To gain a comprehensive understanding of generic QSL 
phases, it is useful to understand the features of exactly solv¬ 
able Hamiltonians which survive the presence of perturbations 
that spoil their exact solvability. Lurthermore, when such per¬ 
turbations are sufficiently strong they can give rise to quan¬ 
tum phase transitions by destabilizing the QSL. Thus these 
systems present microscopic settings to study quantum phase 
transitions out of a QSL phase, an area which is far from well 
understood. Lrom the material perspective, systematic study 
of the influence of such perturbations, which are inevitably 
present in candidate material systems, is also an imperative is¬ 
sue. Motivated by the above questions, in this paper we study 
an example of a concrete spin Hamiltonian that exhibits an 
exactly solvable QSL ground state and additional interactions 
lead to a continuous quantum phase transition to a magneti¬ 
cally ordered state. We systematically study the nature of this 
continuous quantum phase transition out of the QSL. 

The spin Hamiltonian that we study is a variant of the ex¬ 
actly solvable Kitaev model (Lig. 1) on a honeycomb lat¬ 



FIG. 1. (Color online) The Kekule-Kitaev model.^^ The blue (dark), 
green (dashed) and red (light) links are the x,y and z links re¬ 
spectively. There are six sites in the unit-cell as shown (denoted 
by z = 1, 2,3,4, 5,6). We call the sublattices 1,3, 5(2,4,6) as 
odd(even) sublattices. Connecting the mid-point of any particular 
set of links (say links) gives a Kagome lattice (shown in dotted 
fray lines). The lattice vectors are : ai = {3,0}; a 2 = 11, |- 


tice, and we analyze the effect of introducing additional an¬ 
tiferromagnetic Heisenberg interactions and also a magnetic 
fleld. In the present model the distribution of the x, ^ and ^ 
type bonds form a Kekule pattern as shown in Lig. 1. This 
leads to important differences from the original construction 
of Kitaev, with interesting consequences in the structure of the 
phase diagram arising already in the exactly solvable limit. 
In particular, the model reduces to a toric code model on a 
Kagome lattice in an appropriate anisotropic limit. Using a 
combination of analytical and numerical approaches, we show 
that the effect of the Heisenberg interactions are quite differ¬ 
ent in the present case from the by now well known usual 
Heisenberg-Kitaev model. In particular, we describe a mag¬ 
netically ordered phase stabilized by a quantum ‘order by dis¬ 
order’ mechanism^^ and a continuous quantum phase transi¬ 
tion between this ordered phase and a Z 2 QSL in the toric 
code limit. We construct the fleld theory which suggests that 
the critical point belongs to the 3D-XY x Z 2 universality 
class and support this by Monte Carlo simulations. Thus this 
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model presents a controlled microscopic setting for a continu¬ 
ous quantum phase transition between a phase with collinear 
magnetic order and a Z 2 QSL, in itself a subject of much re¬ 
cent interest. 

The rest of the paper is organized as follows. We start with 
a brief introduction to the Kekule-Kitaev model in Sec. II and 
outline the basic features of the exact solution, the phase dia¬ 
gram, and point out the important differences with the origi¬ 
nal Kitaev model. We then discuss an interesting limit of the 
model, the so called “strong-bond limit”, which leads to the a 
toric code model on a Kagome lattice. In Sec. Ill, we add an 
antiferromagnetic Heisenberg interaction to the above Kitaev 
model and investigate the stability of the QSL. While the QSL 
is stable to the weak short ranged spin-spin interactions, as 
expected, they lead to interesting phase transitions when they 
become sufficiently strong. In particular we find that in the 
toric code limit such interactions lead to a magnetic ordering 
that is stabilized by a quantum ‘order by disorder’ mechanism. 
The transition between the QSL and the magnetically ordered 
phase is continuous. Using a combination of field theoretic ar¬ 
guments and Monte-Carlo calculations we find that this con¬ 
tinuous transition belongs to the 3D-XY x Z 2 universality 
class. The effect of an external magnetic field is studied in 
Sec. IV where we find that unlike the usual Kitaev model, the 
present one does not harbour a chiral spin liquid at small mag¬ 
netic field. We summarize our results in Sec. V. Calculational 
details are discussed in the appendices. 


II. THE MODEL 

We start by outlining the spin-1/2 Kekule-Kitaev model on 
the honeycomb lattice. Kitaev, in his pioneering work,^ con¬ 
sidered a spin model on a honeycomb lattice where, depend¬ 
ing on the direction of the three nearest neighbours, there are 
three types of spin exchanges. As pointed out by Kamfor 
et. Kitaev’s original construction of the exactly solv¬ 

able model can be extended to other types of distributions of 
the bond types on the honeycomb lattice. The general Kitaev 
Hamiltonian is given by 

{ij)—a links 

where af {a = x, z) are the Pauli matrices representing the 
spin-1/2 at the site i, and the summation runs over the links of 
the honeycomb lattice, which are of three types (a = x,y, z). 
Here, following Kamfor et. al. we consider a different dis¬ 
tribution of the three types of bonds compared to Kitaev’s 
original model.^ This is depicted in Fig. 1. There are three 
distinct types of hexagonal plaquette, which we denote as: (1) 
x-plaquettes where the bonds alternate between y and ^ types, 
(2) ^-plaquettes where the bonds alternate between x and 2 ; 
types, and (3) 2 ;-plaquettes where the bonds alternate between 
X and y types. The links of a given type are therefore not 
parallel, but instead form a Kekule type of pattern, and so we 
refer to the model as the Kekule-Kitaev model. 

The distribution of links requires the unit cell to contain six 


sites (see Fig.l). We choose the ^-plaquette as the unit cell. 
These plaquettes form a triangular lattice. The Brillouin zone 
information, along with its connection to the Brillouin zone 
of the underlying honeycomb lattice, are given in Appendix 
A (Fig. 9). The symmetries of the above Hamiltonian (see 
Fig. 1) include - (1) lattice translation along the lattice vec¬ 
tors: ai, a 2 , (2) 27r/3 rotation about the plaquette centre, (3) 
refiection about a line connecting the bond centres that lie on 
the opposite side of the plaquette, and (4) time reversal. In 
addition, the model has an extra symmetry along the isotropic 
line Jx = Jy = Jz- This is composed of a simultaneous in¬ 
version of the lattice about an a-bond and a global rotation of 
the spins by tt about the a-axis {a = x^y^ z). 

The exact solution : The exact solution of the above 
model is analogous to that of the usual Kitaev model. ^ Here 
we outline the essential features, and relegate further details 
to Appendix B. As in Kitaev’s original construction we first 
identify conserved plaquette operators. Since the lattice of the 
present model contains three different types of plaquettes, we 
define three types of plaquette operators 

Wc«(p)=-n^?= n 

iG-P ijfG/3—Zin/c,GP 

( 2 ) 

with a = x,y,z. These differ from those introduced in Ref. 
22 by an overall minus sign. By construction, these plaque¬ 
tte operators commute with the Hamiltonian as well as among 
themselves. Hence the Hamiltonian has an infinite set of con¬ 
served quantities which are the Z 2 fiuxes through the plaque¬ 
ttes. These conserved fluxes give rise to fiux sectors, each of 
which has dimension In Appendix B, we show that 

Lieb’s theorem^^ can be used to establish that the ground state 
lies in the zero-fiux sector, where = +1 for all plaque¬ 
ttes. The remaining details of the solution proceed exactly as 
in Kitaev’s construction^ and are outlined in Appendix B. 

The structure of the phase diagram in the present case is 
however quite different from that of Kitaev’s original con¬ 
struction, due to the difference between the symmetries of the 
two models. The phase diagram is shown in Fig. 2, where the 
brown and yellow regions denote gapless and gapped phases 
respectively. The gapless region includes the plane given by 
the equation Jx ^ Jy ^ Jz = In terms of the Majo- 
rana c fermions (see Appendix B), representative dispersion 
curves are presented in Fig. 11. A noteworthy feature is that 
the Majorana c fermions are gapless along the isotropic line 
Jx = Jy = Jz^ as seen in Fig. 2. Along this line, the 
Majorana c fermions of the zero-fiux sector have a nearest 
neighbour tight-binding Hamiltonian with a Dirac point oc¬ 
curring at the F point of the folded Brillouin zone (refer to Fig. 
9), which is is protected by the special inversion symmetry 
present along this isotropic line. Once we move even infinites¬ 
imally away from this isotropic line the Majorana fermions 
gain a mass. The anisotropy so generated is similar to the 
Kekule superconducting order parameter discussed in context 
of graphene.^^ Details on the structure of the mass term for 
the low energy theory are given in Appendix B. 
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FIG. 2. (Color online) The phase diagram of the Kekule-Kitaev 
model (Eq. 1) as a function of the coupling constants. The brown and 
yellow regions represent gapless and gapped phases respectively. For 
isotropic couplings {Jx = Jy — Jz), the system is always gapless, 
but it becomes gapped for infinitesimal perturbations away from this 
limit (see text for details). The band structures for the cut (a-e) are 
shown in Fig. 11. 


A. The strong bond limit: toric code on the Kagome lattice 


An interesting limit of the present model is obtained when 
one of the couplings (say Jz) is much stronger than the other 
two. This is the so called toric code limit.^’^^ Here the Z 2 
fluxes (Wp) provide the low energy degrees of freedom, as 
the c-Majorana fermions (see Appendix B) have a large gap 
(0(Jz)) in this limit.^’^^ To obtain an effective description, 
we first consider the extreme limit Jx = Jy = 0, and Jz > 
0. The lattice separates into disjoint 2 ;-bonds with an Ising 
coupling term —Jzcr^aj. Each bond has doubly degenerate 
ground states | tt)^ I well as two high energy states 

I ti). I it)- We introduce the bond doublets^ 

11) = I tt), U) = I U) (3) 

to represent the ground state subspace on each bond. We adopt 
the notation that the first (second) spin always belongs to a site 
of the sub-lattice of the honeycomb lattice that is denoted by 
open (solid) circles in Fig. 1. It is useful to note that under 
time-reversal symmetry: | \ J|),| J|) ^ | ^). Intro¬ 
ducing Pauli matrices ^ r^) on this bond-doublet space, 

time reversal is affected by T = where X is the com¬ 
plex conjugation operator. This acts as (T : ^ 

—T^}, and so r represent a non-Kramers doublet. 

The toric code appears once small Jx , Jy couplings are 


taken into account, and can be obtained using degenerate per¬ 
turbation theory on the bond doublets. This gives rise to an 
effective model on the Kagome lattice, which is formed by 
joining the midpoints of the z-bonds (refer to Fig. 1). The de¬ 
tails of the effective Hamiltonian are given in Appendix 
giving (up to constants, and at leading nonzero order, up to 6th 
order in perturbation theory): 


Htc — — 


8J| 


+ 


sjl 

8J? 




64J| 


E ^a3-v 

(A,V> 


3^ 

256J| 




( 4 ) 


Here the sum (,) is taken over corner-sharing pairs of trian¬ 
gles, and we have introduced the notations 

= ]^ = ri = ]^ : (5) 

KeA Kev KeO 

where the subscript K refers to the sites of the Kagome lat¬ 
tice. The plaquette operators mutually com¬ 

mute and so the Hamiltonian can be diagonalized simultane¬ 
ously with them. The eigenvalues of the plaquette operators, 
±1, are good quantum numbers, and they specify the ground 
states as well as the excited states. We note that the triangle 
terms do not break time reversal symmetry, although they are 
of the form TjTjT^, because r is a non-Kramers doublet as 
discussed above. While the properties of the toric code model 
are known in context of the square lattice, we shall briefly 
summarize the details here for the sake of continuity. 

The Kagome lattice has a unit cell comprised of three sites 
(for example the up-triangle in Fig. I), and lattice vectors ai, 
a 2 . If Ni and N 2 are the linear dimensions along the two 
crystalline directions, then the total number of up-triangles is 
Nt = A^iA^ 2 - There are 3Nt Kagome lattice sites, and thus 
the Hilbert space spanned by the r spins has a dimension of 
23ArT Correspondingly there are Nt of each of the Ta, Ty 
and Xq operators, and specifying their values also leads to 
2^^^ states. However, the operators are not all independent 
on a 2-tori, as they obey two constraints: 

(nM(nA)=+i. n 9^0 = +1 • (6) 

V A O 

These give rise to the topological degeneracy of 4, as expected 
for a gapped Z 2 quantum spin liquid in a toric code model. 

On taking Jy = 0, the lattice splits into disconnected up- 
triangles. There are eight states per triangle with a four-fold 
ground state degeneracy. The ground state degeneracy for a 
lattice of 3Nt sites is then 4^^ = 2^^^. This extra degener¬ 
acy is accidental and is immediately lifted when Jy is turned 
on. For < 0 (we shall always take Jx > 0), the ground 
state lies in a sector where the eigenvalues of each of the pla¬ 
quette operators are Ta = X^ = —Xq = +1, as this min¬ 
imises each of the four terms in the Hamiltonian. For even 
Nt, the ground state has the 4-fold topological degeneracy 
referred to above. For odd Nt however, the constraints of 
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FIG. 3. (Color online) (a) The medial honeycomb lattice (red) of the 
Kagome lattice (blue), (b) A boundary can be obtained by setting the 
coupling strengths along the red (dashed) triangles and honeycombs 
to zero. 


Eq. (6) do not allow all hexagons to have = —1, and 
so the ground state must have one of them to be +1. This 
“defect" honeycomb can sit anywhere on the 2-tori and the 
ground state degeneracy is raised to ANt. A similar feature is 
also observed for > 0. 

In the bulk there are two types of excitations, both of 
which are gapped and dispersionless. ^ There are Ising elec¬ 
tric charges which are associated with the triangular plaque- 
ttes, and Ising magnetic charges associated with the hexagonal 
plaquettes. To study these excitations it is instructive to go to 
the medial lattice of Kagome, which is obtained by joining the 
centres of neighbouring triangular plaquettes as shown in Fig. 
3(a). This gives another honeycomb lattice, which is different 
to that of Fig. 1. The excitations can now be understood as 
the endpoints of string operators. A pair of electric charges 
are created by the operator 

(7) 

lei 

where i denotes a path on the honeycomb lattice starting and 
ending at the two charges, while a pair of magnetic charges 
are created by 

( 8 ) 

lei 

where here ^ is a path on the triangular lattice obtained by con¬ 
necting the centres of the honeycomb plaquettes. Examining 
how the charges wind around one another, it is straightforward 
to see that the electric and magnetic charges are bosons with 
mutual semionic statistics, i.e. each sees the other as source 
of TT flux.^ We note that while the electric charges move on a 
bipartite (honeycomb) lattice, the magnetic charges move on 
the non-bipartite (triangular) lattice. 

Open boundary conditions and Majorana edge modes: 
The system has gapless edge states in the presence of an 
open boundary, which are related to the underlying topologi¬ 
cal order. To describe these, let us imagine drawing a great 
circle around the system on a 2-tori, obtained by translation 
along ai as shown in Fig. 3(b), and then cutting the sys¬ 
tem by setting all the ITa , crossed by the great 

circle to be zero. Let us again take and N 2 as the lin¬ 


ear dimensions, so there are 3A^i A ^2 — connected Kagome 
lattice sites, and the numbers of each of up and down trian¬ 
gles and hexagons are A^i(A ^2 ~ 1)^ leading to a degeneracy 
of On the other hand there 

are 4A^i edge spins, 2A^i on each edge. This gives rise to 
a/2 degree of freedom at each edge site, which corresponds 
to Majorana edge states. These edge states have a flat band 
with exactly zero energy and are stable to weak perturbations 
away from the exactly solvable limit as long as translational 
symmetry along the edge is not spontaneously broken. There 
are other more complicated edges that could be considered, 
but this is outside the scope of the present work. 

Distinct toric codes and adiabatic continuity: In the 
above treatment we have focused on the toric code model that 
appears in the Jz /c, Jy regime. However we could have 
done the same for the other two couplings, Jx and Jy. In Ki- 
taev’s original model^ the three toric codes so obtained cannot 
be connected without closing the bulk excitation gap. How¬ 
ever, we And that in the present model there exist only two 
such distinct toric code limits. These are separated by the 
Jx A- Jy A- Jz =0 plane and are adiabatically connected to 
the two toric codes obtained in the vicinity of the two limits 

Jz —^ dz(X). 

To demonstrate this, imagine setting Jx = Jy = 0 and 
Jz = I, which breaks the honeycomb lattice of Fig. 1 into 
disjoined bonds (which are the sites of the Kagome lattice). 
Turning on small Jx only results in disjoined honeycombs 
(or isolated up-triangles for the Kagome lattice). This is 
contrary to the original Kitaev model where turning on two 
bonds results in extended ID chains. Thus in the present case, 
even for flnite Jx, there is no dispersion when = 0. For 
the disjoined honeycombs, albeit with anisotropic exchanges 
on their bonds, the six eigenstates have energies Az{Jx + 

J,), ±v/J| -J,Jz + j|, ±^/Jl - j.j. + . For 

any value of > 0 the gap never closes. The wavefunctions 
however become superpositions of the a spin states on the 
six sites of the hexagon containing x and 2 ; bonds. On 
attaining the point Jx = Jz = one can then take Jz to zero 
gradually without closing the gap, resulting again in disjoined 
bonds, but now the doublets are polarized parallel to the x 
direction. So it is possible to continuously connect the two 
bond limits without closing the bulk energy gap, and hence 
to adiabatically connect the two toric codes which arise from 
perturbations about these two limits. This argument does not 
work however if we try and connect Jx = Jy = 0, Jz = ^ 
io Jz = Jy = Jx = —1, as then a gap closing point is 
necessarily encountered. It follows then that the two toric 
code phases obtained in the vicinity of > 0 and Jz < ^ 
cannot be continuously connected as the bulk gap closing and 
reopening necessarily occurs across the Jx + Jy + Jz = 0 
plane. 

This completes our discussion of the phase diagram of the 
Kekule-Kitaev model. Starting from the next section we shall 
investigate the effect of perturbations. 
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III. HEISENBERG PERTURBATIONS 

An important class of interactions in effective spin models 
are bilinear spin-spin interactions. To this end we now inves¬ 
tigate the effect of adding antiferromagnetic Heisenberg in¬ 
teractions to the Kekule-Kitaev model, and shall see that this 
leads to interesting consequences. We thus turn our attention 
to the following extended Hamiltonian 

H = Jh ^ (Ji ' aj ^ IKk , (9) 

(*i> 

where IKk is given by Eq. 1 . As we shall see, the resulting 
phase diagram is quite rich and allows for interesting phases 
and phase transitions, the most interesting of which is a con¬ 
tinuous transition between the QSL and a magnetically or¬ 
dered phase. In this paper, we shall concentrate on the toric 
code limit for the Kitaev interactions and obtain a controlled 
description of not only the magnetic ordering driven by the 
Heisenberg perturbation, but also of a continuous quantum 
phase transition between the QSL and the magnetically or¬ 
dered phase. 


A. The generalized toric code model 

Starting with Eq. 9, in the limit Jn^Jx^Jy ^ Jz^ we ob¬ 
tain a generalized toric code model with nearest neighbour 
antiferromagnetic Ising perturbations. The hierarchy of en¬ 
ergy scales allow for a strong coupling perturbative expan¬ 
sion, as in the previous section, yielding the following effec¬ 
tive Hamiltonian to the leading order in all the couplings 

H =a'^ r|rj - oa E] S'a - Ov X] “ «0 X] 9“o 
(IJ) A V O 

-flAv 51 3'A^fv (10) 

(A,V> 

This is a generalized toric code model where the last four 
terms are same as the toric code Hamiltonian (Eq. 4), albeit 
with renormalized couplings, and the first term is a nearest 
neighbour antiferromagnetic Ising interaction between the r 
spins sitting on the Kagome lattice, and a = Jh to leading 
order. 

A special feature here is that remains a conserved quan¬ 
tity, which is not the case in the original version of the Kitaev 
model.As a result, the magnetic fluxes are still good quan¬ 
tum numbers. This will be very important for the rest of our 
analysis and will help reveal the nature of the phase transi¬ 
tion in the present model, which has so far eluded analytic 
understanding in the original Kitaev-Heisenberg models even 
in the toric code limit.The perturbations do however cause 
the electric charges to acquire dynamics. 

To see the effect of the Ising term more clearly, it is useful 
to turn once more to the medial lattice construction shown in 
Eig. 3(a). The resulting honeycomb lattice (not to be confused 
with the lattice of Eig. 1) has two sublattices which are respec¬ 


tively at the centres of the up and the down triangular plaque- 
ttes of the Kagome lattice. The Ising term leads to hopping of 
electric charges on the honeycomb lattice preserving the sub¬ 
lattice fiavour, i.e. next-nearest-neighbour hopping. Since the 
electric charges are bosons, they can condense once their dis¬ 
persion minimum touches zero. We shall show that this phase 
breaks time reversal symmetry and hence generates magnetic 
order. 

Let us first generalise our viewpoint. While the Hamilto¬ 
nian of Eq. (10) has been derived from a microscopic model 
with a definite hierarchy of the coupling constants, we shall 
relax that hierarchy for the moment and consider the general¬ 
ized phase diagram for the above model. We shall comment 
on the actual microscopic parameters at the end in context of 
the generalized phase diagram. 

To proceed we introduce a gauge theory description of 
the generalized model, by defining the following set of Ising 
variables^ ^ on the medial honeycomb lattice (Eig. 3(a)) 

3“a=Mp, (11) 

where /i spins are defined on the sites and are defined on the 
links. The = ±1 carry Ising gauge charge and = ±1 are 
the Ising gauge potentials. Here p, q are the nearest neighbour 
sites on the medial honeycomb lattice, and so belong to the 
two different sublattices of it. The Hamiltonian of Eq. 10 now 
takes the form 

H = a - flA 53 “ “v 53 1^9 

{{pp')),q P q 

- ao 53 n Pp‘i~ 53 

O pqeO (pq) 

where in the first term q is the common nearest neighbour site 
connecting the two second nearest neighbour sites The 

Ising magnetic fluxes are given by 

= (13) 

pqeO 

and they sit on a triangular lattice formed from the plaquettes 
of the medial honeycomb lattice. 

In the limit ua, ^ and a, uq > 0, the Hamil¬ 

tonian becomes classical as it has only operators. The 
most degenerate point of parameter space is obtained on fur¬ 
ther setting ao = 0. This corresponds to a classical anti¬ 
ferromagnetic Ising model on the Kagome lattice, and has 
a ground state entropy of 0.502 ks per site of the Kagome 
lattice.The degeneracy is partially lifted however by an 
infinitesimally positive ao- This chooses the ground state 
sector which has zero magnetic flux through the hexagonal 
plaquettes. This is easiest to see in the gauge theory, where 
ao forces the magnetic flux through the hexagons to be zero. 
As a result, a gauge with all = +1 can be chosen, yield¬ 
ing two copies of the triangular lattice antiferromagnet, each 
of which has an entropy of 0.323 ks per site of the trian¬ 
gular lattice.As there are three Kagome lattice sites for 
each triangular lattice site, the resulting entropy is 0.215 ks 
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per site of the Kagome lattice. The quenching of the entropy 
comes from the fact that spin configurations whose product 
around the hexagons is — 1 are now energetically more costly 
and hence not in the ground state manifold. The entropy nev¬ 
ertheless remains macroscopic due to the absence of quantum 
fluctuations, to which we now turn our attention. 

To consider the quantum fluctuations, let us start with the 
limit ua, ^ ^ with all being positive. In this 

limit we can concentrate on the sector where there are no mag¬ 
netic charges, i.e. Tq = -f 1 Vo. Hence we can choose a 
gauge where = 1 for all the links of the medial lattice. The 
Hamiltonian in Eq. (12) can then be cast into a transverse held 
Ising model (TFIM) 

HtFIM =CL ^ PpPp/ - QA 

{{pp')) p 

~ • ( 14 ) 

{pq) 

In terms of the new variables, the spin liquid of Sec. IIA 
becomes the paramagnetic phase of the p spins in which all 
the spins are polarized along When a = uav = 0 
is clearly the ground state, and the basic excitations, the Z 2 
electric charges (the flipped spins), are gapped with an en¬ 
ergy cost of the order a a or . When a is turned on these Z 2 

charges acquire dispersion, with a bandwidth proportional to 
a. Once this becomes comparable to the gap, these charges 
condense at a particular momentum. In such a condensed 
phase (p^) 7 ^ 0, and consequently (r^) 7 ^ 0. Since is odd 
under time reversal, this phase breaks time reversal symmetry 
(along with lattice symmetries) and is actually a magnetically 
ordered state. 

The magnetic order: Let us first describe this magnetic or¬ 
dering for uav =0’ this decouples the Hamiltonian to two 
copies of the TFIM on the triangular lattice, a model which is 
well known.In the classical limit this model has an ex¬ 
tensive ground state degeneracy, albeit with power-law spin 
correlations.^^ When the transverse held is turned on, the sys¬ 
tem goes into a magnetically ordered state through a quantum 
‘order by disorder’ mechanism.^^’^^ The ground state has a 
characteristic x ^/3 order where there is a hexagonal lat¬ 
tice superimposed on the triangular lattice on which the spins 
exhibit Neel order in direction, while the spins at the site in 
the centre of each hexagon are polarized along the p^ direc¬ 
tion, and gain their energy from the transverse held by being in 
the maximally flippable state. For each triangle the three sites 
have magnetization, Mz = (p^), of the form (+ 1 ,— 1 , 0 ). 
This candidate magnetic order in terms of the p^ spins for a 
single triangular lattice is shown in Fig. 4(a). The ordering in 
terms of the spins can be obtained through Eq. 11 and this 
is shown in Fig. 4(b). A characteristic feature of the order¬ 
ing is the regular pattern of bow-ties with zero magnetization, 
which are surrounded by zig-zag chains of antiferromagnetic 
order. It is interesting to note that the chains of antiferromag¬ 
netic order are mutually decoupled from each other by the 
zero magnetization bow ties. In our Monte Carlo studies on 
the Hamiltonian 14 (described below) on a space-time lattice, 
we And the above ordering pattern persists throughout the part 


of the phase diagram that is magnetically ordered even when 
the triangular lattices are coupled (uav 7 ^ 0 )- 

Now we describe the transition between the QSL and this 
magnetically ordered state, which as we remarked earlier 
should be looked upon as a transition arising from the con¬ 
densation of the Ising electric charge excitations of the QSL. 
Considering again one of the triangular sublattices, it is well 
known^^ that the transition of the TFIM can be effectively de¬ 
scribed by adopting a soft-spin description for the p^ spins, 
and identifying the soft modes that condense to give rise to 
the magnetic order. We follow a similar prescription for de¬ 
scribing the present phase transition. The soft modes occur at 
Ka = ±[ 47 r/ 3 , 0], and so the spins expand as 

p^(r) = V^+(r)e^^+-^ -h V^_(r)e^^--^ , (15) 


with amplitudes (?/;+,?/;_). Taking all the global symme¬ 
tries of the microscopic model into account, the (2 -f- 1 )- 
dimensional Euclidean Landau-Ginzburg action can be con¬ 
structed in the standard way (details are relegated to Appendix 
D) giving^^ 


= J -h r2'0 •-0 + •'0)^ + 

+'^6('0++'0-)] , (16) 


where the integration is over the (2 -f 1)-dimensional Eu¬ 
clidean space-time, and 


7^= 


1 

2 ’ 2i J 


(17) 


The above action predicts^^ that the critical point belongs to 
the 3D-XY universality class, with the six-fold anisotropy 
term being dangerously irrelevant at the critical point. The 
sign of this anisotropy term determines the nature of the mag¬ 
netic order. 


For a A = Uy, the Hamiltonian has an inversion symmetry 
that exchanges the two sublattices of the medial honeycomb 
lattice. In the rest of our calculations, we shall restrict our¬ 
selves to this inversion symmetric case and use = 

r. The limit uav =0 corresponds to two such decoupled 
triangular lattices. For each copy we introduce a pair of soft 
modes given by where n = 1 , 2 denotes the two copies, 
and so the action is 

2 

^decoupled (18) 

n=l 

The six-fold anisotropy term is again dangerously irrelevant 
at the critical point, and the symmetry is 3D-XY x 3D-XY. 
In this limit the magnetic orders of the two triangular lattices 
are mutually independent. Symmetry however allows the cou¬ 
pling of the two triangular lattices and this is exemplifled by 
the presence of the term uav 7 ^ 6 the generalized toric 
code model. We can write down the leading symmetry al¬ 
lowed term (see again Appendix D for details) to get the com- 
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(a) (b) 

FIG. 4. (Color online) (a) Magnetic ordering for the antiferromagnetic TFIM on the triangular lattice. The dots refer to sites which have zero 
magnetization (the maximally flippable sites) (b) Magnetic ordering of the spins on the Kagome lattice. 


plete action 


S — ^decoupled ^int 5 ( 19 ) 

where to the leading order 


Sint — 


J (?’2 


V4 

n=l,2 


( 20 ) 


The term Sint breaks the 3D-XY x 3D-XY symmetry down 
to 3D-XY X Z 2 with a six-fold anisotropy term. The Z 2 sym¬ 
metry is related to the fact that under inversion about the bond 
centre, the flavours of the soft-mode change. Tracing back to 
our microscopic Hamiltonian, this symmetry is present when 
a A = tty- Power counting, at the free flxed point, shows 
^2 ? ^4 5 ^4 ^4 relevant while the sixth order terms {uq 

and vq) are marginal. To understand their effect for the or¬ 
dered phase (r 2 < 0) we look at the phase fluctuations. To 
this end we write 


= p(cos sin^^), (21) 


where we have taken the magnitude of the XY order parameter 
to be constant. Then, neglecting the six-fold anisotropy term, 
for the ordered phase we have^^ 

Sdecoupled = / dv^^ + \de.\^] , (22) 

where 6± = {6i ± 02 )/ 2 , and 


Sint = / d-^T 


I' 


r'n cos 20_ + — cos 46>_ 


(23) 


Now, depending on the sign of r'^ and w^, an expectation 
value for is chosen. Under inversion -> —0-, and so 
a transition which induces an ordering in 0- simultaneously 
breaks the Z 2 symmetry related to inversion. Expanding the 


cosine terms around the 0- expectation value makes 0- ex¬ 
citations massive while remain massless. The latter is in 
fact the low energy XY degree of freedom. Now the six fold 
anisotropy term becomes 

Vq [cos[ 6(6>+ -h 6>_)] + cos[6(6>+ - 6>_)]] 

= 2vq cos(66>+) cos(66>_) . (24) 

Replacing with its expectation value, we And that this 
term behaves like the six fold anisotropy fleld for 6>+, which 
is known to be relevant in the low temperature ordered 
phase^^’'^^’'^^ . Hence we expect that the six-fold anisotropy is 
dangerously irrelevant at this transition, and that the critical 
point belongs to the 3D-XY x Z 2 universality class. This ex¬ 
pectation is supported by our numerical studies below. Thus 
the fleld theory predicts that the transition due to the con¬ 
densation of the Ising electric charges in the QSL belongs to 
an interesting universality class, and that the present micro¬ 
scopic model can harbour such an unconventional phase tran¬ 
sition from a topological QSL phase to a magnetically ordered 
phase. 

Numerical calculations: To complement the above pre¬ 
diction of the soft mode analysis, we investigate the lattice 
model numerically. In particular we construct a discrete-time 
classical action for the spin model of Eq. (14), on which we 
perform Monte Carlo simulations. The details of the deriva¬ 
tion are given in Appendix E, and the resulting three dimen¬ 
sional Euclidean Landau-Ginzburg action, E, describes an 
Ising lattice model on a stacked honeycomb lattice with ad¬ 
ditional Ising link flelds. This has the form: 


JE=J Y 

({ii».T h-r \ j J (25) 

-^T ^ ^ Vij,T * 


where ±1) and r]ij{= ±1) are Ising variables sitting 

on the sites and links of the stacked honeycomb lattice re¬ 
spectively. The couplings are related to those of the micro- 
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FIG. 5. (Color online) The phase diagram of the model (25) in the 
1 / J-KrjJ plane. The phase boundary is obtained from the peaks of 
the specific heat (blue dots), see Fig. 6(a). To ensure that these peaks 
correspond to the phase transition, critical couplings extracted from 
crossings of the Binder ratio are plotted (red crosses) for Kr j J — 
2,4,6, see Fig. 7(a) for K^jJ — 2. For small K^jJ the peak 
in the specific heat becomes rounded and the accuracy of the phase 
boundary is diminished. The two dashed lines represent the curves 
along which the magnetic correlations are plotted in Fig. 6(b). 


scopic model as follows: J = aAr, tanh(iC) = , 

tanh(AraAv) = (thus K^- represents the coupling 

between the two triangular sublattices), and Ar = P/N, 
where p is the inverse temperature (we are interested in the 
P ^ oo limit) and N is the number of temporal slices. 
Without loss of generality in the scaling limit, we study the 
isotropic case iC = J for the above model and explore the 
phase diagram as a function of the two dimensionless param¬ 
eters 1/J and KrjJ. The coupling J controls the strength of 
the Heisenberg perturbations, while controls the coupling 
of the two copies of TFIM. In addition, 1/J plays the role of 
temperature for the classical action. 

In Fig. 5, we plot the phase diagram as a function of 1/ J vs. 
K^l J, which is obtained from the peaks of specific heat (Fig. 
6 (a)), and complemented by the calculation of selected criti¬ 
cal couplings from the crossings of Binder ratios. We remind 
the reader that the paramagnetic phase of above model of Eq. 
(14) corresponds to the spin liquid phase of the microscopic 
Hamiltonian. 

To explore the stability of the spin liquid to small Heisen¬ 
berg perturbations we plot the specific heat against 1/J for a 
range of values of J in Fig. 6 (a). This shows that the 
spin liquid is stable for small J, and that it remains so until 
a critical value is reached, with the sharp peak in the specific 
heat (and concomitant development of magnetic correlations) 
indicating a phase transition to the magnetically ordered state. 
We remark that the broad feature of specific heat seen at large 
1/J for large Kj- in Fig. 6 (a), which is absent in the Kj- oc 
limit^^, is due to the thermal excitation of the link variables 
r]ij^r which occurs at 1/J ^ K^j J. 

Below the critical values of 1/ J the system is magnetically 








1/J 



KJJ 


(b) 


FIG. 6. (Color online) (a) Plots of the specific heat c as a function 
of 1/ J for a range of values of KrjJ. (b) The sublattice magnetic 
correlations plotted along the lines KrjJ — 2 and 1/J = 2. Both 
sublattices have an identical ordering as described in the text, and 
there are no magnetic correlations between the sublattices. 


ordered as described before. To investigate the nature of the 
ordering, we take a 3-site unit cell for each sublattice, and in¬ 
vestigate the correlations within a sublattice, and between the 
two sublattices. The correlations within a sublattice are plot¬ 
ted in Fig. 6 (b) for two cuts of the Ij J-K^ plane, given by 
Kr!J = 2 and IjJ = 2 respectively. For these the mag¬ 
netisations are ordered as Mi > M 2 > M 3 at each step of 
the Monte Carlo simulation, indicating the magnetic ordering 
(M|) - 0, {Ml) - (M|) - -(M1M3) displayed in Fig. 
4(a). For the computations described in this section we take 
the order parameter to be M = Mi — M3. We use this to 
calculate the Binder ratio R 2 = (M^)/(M^)^ to determine 
the position of the phase transition to higher accuracy. This 
is shown in Fig. 7(a). We find no correlations between the 
respective orderings of the two sublattices, yielding the full 
magnetic ordering of Fig. 4(b) in terms of the r spins. From 
Fig. 6 (b) it is clear that the magnetic order does not change 
within the ordered region as far as we can resolve within our 
numerical calculations. 

To probe the critical point further, we study the critical ex¬ 
ponents, which we extract via a scaling collapse of the mag¬ 
netic susceptibility as shown in Fig. 7(b). We restricted 
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1.97 2 2.03 2.06 2.09 2.12 


1/J 

(a) 



(b) 


FIG. 7. (Color online) (a) The specific heat c near the transition for 
KrjJ — 2 and the Binder ratio gives the value of the critical 
coupling Xj Jc — 2.0505. (b) The scaling collapse of the magnetic 
susceptibility. 


to small lattices, and focused on L x L x L with L = 
9,12,15,18, 21, 24. These rather small sizes to which we are 
restricted lead to notable finite size effects. Within error how¬ 
ever, we find that the scaling of the critical exponents for finite 
Kr match those for = oo. In addition, for = oc we 
went to L = 27, 30 and obtained results consistent with the 
3D-XY critical exponents. While we do not completely un¬ 
derstand why the present exponents are very close to those for 
3D-XY, to the best of our knowledge the actual exponents for 
a 3D-XY X Z 2 critical point are not known. A related issue 
is whether one can separate the 3D-XY and the Z 2 transi¬ 
tion to open up a phase in between where the Z 2 (inversion) 
symmetry is broken but where there is no magnetic order. We 
have not been able to achieve this in the present model, but 
this may be due to the fact that the Z 2 can be described by 
Ising variables sitting on the links of the medial honeycomb 
lattice, which is the Kagome lattice. This may cause frus¬ 
tration among the Z 2 variables and could prevent them from 
ordering alone in absence of the XY field. 


The parameters of the generalized toric code and the mi¬ 
croscopic model: we conclude with some comments on the 
relationship of the parameters of the generalized toric code 
Hamiltonian (Eq. (10)) to the parameters of the microscopic 
model of Eq. (4). Eirstly, the couplings ua and neces¬ 
sarily appear with a relative negative sign in Eq. (4) when 
oq > 0, whereas we consider a a = relative sign 

can be removed however by a transformation which rotates the 
p spins about the axis on one of the triangular sublattices 
of the medial lattice such that for that sublattice. 

Next, we comment that the presence of the Heisenberg inter¬ 
action would cause the coupling constants to change. These 
renormalizations would be small when Jx^Jy > Jh (we as¬ 
sume Jx^Jy > 0). Einally, we have done our calculation in 
the zero magnetic charge sector which is the relevant sector in 
the regime where uq is the dominating coupling constant of 
the generalized toric code model (eq. 10). This may not be so 
however as suggested from the couplings of the microscopic 
model. In the present case we assume that of the different 
magnetic charge sectors (which remain good quantum num¬ 
bers in presence of the Heisenberg perturbations), the zero 
charge sector always remains the ground state. Numerical ver¬ 
ification of this assumption forms a topic of future study. 


IV. EFFECT OF A MAGNETIC FIELD ON THE 
KEKULE-KITAEV MODEL 

We now briefly discuss the effect of a magnetic field on the 
Kekule-Kitaev model. We consider a Hamiltonian of the form 

Hz = Hk - (fi (26) 

where h is the magnetic field. As a single spin operator creates 
two fluxes which cost energy, in the limit of small fleld the 
effect of the time-reversal symmetry breaking Zeeman term 
can be obtained by perturbation theory within the ground state 
(zero flux) sector. 

Let us thus describe the effect of the Zeeman term on 
the Majorana fermions. The first non-trivial interaction that 
breaks time reversal symmetry is a three spin term, which 
when written in the zero flux sector provides next nearest 
neighbour hopping to the Majorana fermions (similar to the 
case in the original Kitaev model, other terms renormalize 
the nearest neighbour hopping or provide short range four 
fermion interactions which are irrelevant at the free Majorana 
fixed point). At the isotropic line (Jx = Jy = can use 

a two site unit cell (used only in this section) for the Majorana 
fermions to make our results transparent. The hopping Hamil¬ 
tonian obtained from including the Zeeman term through per¬ 
turbation theory is shown in Eig. 8. It must be noticed that 
as one goes along the arrows in one hexagonal plaquette, the 
winding of the red and the blue arrows are mutually opposite 
contrary to the case in the usual Kitaev model.^ This means 
that the mass term obtained at the Dirac point does not invert 
from one valley to another and hence the Chern number of the 
two gapped bands are zero. This means that the state so ob- 
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avenue of future study. 
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Appendix A: The details of the lattice and Brillouin zone 


FIG. 8. (Color online) The nearest and the next nearest neighbour 
hopping for the isotropic model with two site unit cell. Going along 
the arrow, connoting two sites i and j gives a term iacj in the tight 
binding Hamiltonian. 


tained is not a chiral spin liquid as was obtained in the original 
Kitaev model^ and hence does not support gapless edge modes 
in this gapped time reversal symmetry broken QSL state. 


V. SUMMARY AND DISCUSSIONS 


Motivated by the idea to understand quantum spin liquid 
phases, and phase transitions from them, we have studied a 
variant of the Kitaev model with and without antiferromag¬ 
netic Heisenberg interactions. We found that the phase di¬ 
agram even in the absence of perturbations is quite different 
from the original Kitaev model, and includes regimes captured 
by toric code models on the Kagome lattice. In the presence 
of Heisenberg interactions, in the toric code limit, the system 
shows an interesting quantum phase transition to a magneti¬ 
cally ordered phase where the magnetic order itself is chosen 
through a quantum ‘order by disorder’ mechanism. Using a 
combination of field theory and Monte Carlo studies we have 
been able to study the phase transition which likely belongs 
to the 3D-XY x Z 2 universality class. Such a controlled de¬ 
scription of a continuous quantum phase transition out of a 
QSL to a non-trivial magnetically ordered phase in a micro¬ 
scopic model is quite interesting in the context of recent in¬ 
terest in understanding novel QSL phases. We note that it has 
been claimed"^^ that a transition belonging to the 3D-XY x Z 2 
universality class can be driven by fluctuations to a weak first 
order transition. However, we have not seen signatures of such 
a discontinuous transition in our numerics. While the possibil¬ 
ity of a weak first order transition cannot be completely ruled 
out, we hope that future numerics on larger system sizes will 
be able to solve this issue. 

The nature of the phases and phase transitions arising from 
Heisenberg antiferromagnetic perturbations in the limit of 
isotropic Kitaev couplings appears to be more complicated, 
due to the absence of a special SU (2) invariant point^^ in 
the present case. While from general considerations the QSL 
phase is expected to be stable to Heisenberg perturbations, 
a determination of the window of stability and the resultant 
phase to which this QSL gives way constitutes an interesting 


The lattice vectors, as shown in Fig. 1, are given in the 
caption of the same figure. The Brillouin zone is presented in 
Fig. 9, along with the reciprocal lattice vectors. 



FIG. 9. (Color online) The Brillouin zones: The red (black) lines 
and arrows denote the Brillouin zone of the isotropic (anisotropic) 
honeycomb lattice with two (six) point unit cell. The Reciprocal 
vectors are denoted by Pi{Gi) and P 2 {G 2 ) for the two (six) site 

cases. Gi = G 2 = |o,^|- We note that the 

zone corners of the two-site unit cell lattice maps to the zone centre 
of the six-site unit cell case. 


Appendix B: The exact solution of the Kekule-Kitaev model and 
the applicability of Lieb’s theorem 

To obtain the exact solution, following Kitaev, we can de¬ 
fine the spins in terms of four mutually anticommuting Majo- 
rana fermions 


af = ih^Ci (Bl) 

subject to the constraint 

D, = b^b^,b^,Ci = l, Vi. (B2) 

The Di commute with the Hamiltonian, [H^Di] = 0, re¬ 
sulting in a Z 2 gauge structure, for which the Di generate 
the Z 2 gauge transformation. Continuing a similar treatment 
as in the original case we write for the a-th link crfcr" = 
{ibfci){ib'jCj) = -iufjCiCj, where, rtg- = ibfb'j = —u^i- 
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Generation of mass around the isotropic line — Jy — Jz) 


On the isotropic line, we have four bands, with linear dis¬ 
persion, touching at the V point. The low energy k.p Hamil¬ 
tonian for these four bands (for c majorana fermions) can be 
obtained (for the zero flux sector) from the lattice hopping 
Hamiltonian by expanding about the T point to the linear or¬ 
der and projecting to the four bands. This has the form 


^iso = 


3 

4 


0 0 kx 

0 0 e^'^^ky 

p Jc p ^^3 n 

P ^^2 p il~ 4 . Q 

o ^ ^x ^ 


pi'T2 Ip 
o ^y 

pii~4 
o ^x 

0 

0 


(B5) 


FIG. 10. (Color online) The Kekule-Kitaev model on a brick lattice, 
which is a square lattice with selectively deleted bonds. The reflec¬ 
tion symmetry about the black line allows the application of Lieb’s 
theorem. 


All the ufj commute amongst themselves as well as with the 
Hamiltonian, and are thus constants of motion. However they 
are not however gauge invariant. Rather they are the gauge 
potentials of the Z 2 gauge theory. The plaquette operators of 
Eq. (2) take the form 


w„(p) = n 

ij£/3—link,EP 

and the Hamiltonian becomes 


^ a; i G odd^j G even) 

(B3) 

(B4) 




a—links 


where (we have numerically calculated) ri = 
0.8714997r,r2 = -0.13827l7r, T 3 = -0.08333337r, T 4 = 
—0.09310377r. The doubly degenerate eigenvalues are given 
by C±(k) = ±||k|. This is the linearly dispersing band 
touching shown in Fig. 11 . 

We now wish to move away from the isotropic line by 
adding an anisotropy of the following form: = J ^6^ Jy = 

Jz = J. This anisotropy pattern is similar to the Kekule order 
parameter of Ref. 27 (with a = 0 in their notation in Eqs. 9 
and 10 of that paper). The mass matrix has the form: 


^mass 


6_ 

4 


0 0 
0 0 


\/3e*^3 

0 

0 


\/3e*^2 - 

0 

0 


(B6) 


where r]i = 0.3714997r, 7^2 = 0.3617297r, 7^3 = 

—0.5833337r, 774 = 0.4068967r. The corresponding doubly 
degenerate bands have the dispersions. 


We can replace the with their eigenvalues ±1, recasting 
the systems as a tight-binding model where the c Majorana 
femions hop on a bipartite structure such that the hopping am¬ 
plitudes tij 7 ^ 0 only if 7 G odd and j G even, and all 
loops contain an even number of edges and are planar. Keep¬ 
ing intact the relevant symmetries of the model, the lattice can 
be deformed to the one shown in Fig. 10. Clearly, the figure 
shows that the results of the Lieb’s theorem^^ are applicable 
in this case and hence the ground state lies in the zero flux 
sector for the Majorana fermions. This fiux sector is ensured 
by choosing u^- = when i G odd and j G even. Let us 
note an important difference with the usual Kitaev model. As 
each plaquette contains a definite kind of bond three times or 
not at all, a global sign change in any of Jx , Jy or Jz can no 
longer be absorbed in changing the sign of the corresponding 
u^j (through a gauge transformation). As a result the ground 
state energy in the sectors {Jx^Jy^Jz) and {Jxi Jy^—Jz) will 
be different. The band structure of the Majorana fermions 
is obtained by diagonalizing the free Majorana Hamiltonian, 
giving rise to 6 particle-hole symmetric Majorana bands. We 
plot them and analyze the band structure in different parame¬ 
ter regimes to find the phase diagram shown in Fig. 2. Some 
representative band structures are shown in Fig. 11. 


C±(k) = ±y2k2 + i52, (B7) 

The double degeneracy is lifted away from the Dirac point due 
to higher order terms as seen in the Fig. 11 . 


Appendix C: Perturbative derivation of the toric code model 

The toric code Hamiltonian Htc of Eq* (4) appears 
through degenerate perturbation theory on restricting the 
Hilbert space of the full model to the ground state manifold 
of the bond doublets | 1)^), | J|) of Eq. (3) at the z-bonds of 
the honeycomb lattice (which form a Kagome lattice). The 
Hamiltonian is computed order by order in 1/ 

H^^f=P[V^ VG'qV + VG'qVG'qV + • • • ] P (Cl) 

where G'o = (1 - P)Pjrp - P) with P = 

YIk (I (iT I + I ^)kk (J| I) being the projector onto the 
ground state manifold {K labels the sites of the Kagome lat- 
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FIG. 11. (Color online) The Majorana band structure: (a)-(e) denotes the cut shown in Fig. 2. (f)-(j) denotes the low energy section of the 
bands corresponding to (a)-(e) respectively. 


lice), and 

z —links 

V=-Jy Y Y (C3) 

y —links x—links 

The first non-trivial terms occur at third order 


where ^nd (,) indicates the sum is to be 

taken over neighbouring pairs of triangles. The operators Jq 
are equivalent to the MOz plaquette operators of the full model, 
while the operators 9^a and are equivalent to the Wy and 

(3) 

Wx plaquette operators respectively, and so Htc = ^e// + 

is a good effective model, it captures the leading physics 
of all degrees of freedom in the strong bond limit. 


^(d) _ 
- 


o t3 

^ A,{I,J,KeA} 

Q /3 

ot2 a 

8 J? rj2 a 


A 


V ’ 


(C4) 


while next non-trivial contribution to the effective Hamilto¬ 
nian occurs at sixth order. The second and fourth order terms 
are constant, while the fifth order term gives a subleading cor¬ 
rection to the third order term. At sixth order there is a new 
term associated with each hexagon of the 

Kagome lattice, where Ij labels the z-bonds at the vertices 
of the hexagon. In addition there is a contribution from pairs 
of neighbouring A, v triangles, where 

Jj , Kj label the z-bonds at the vertices of the two respective 
triangles, one of the bonds being common to both. Counting 
up these two contributions gives the relevant sixth order term 
of the effective Hamiltonian 


Me) _ 

^eff - 


3/3/3 

256/3 


o 


7/3/3 


64/3 


(A,v> 


(C5) 


Appendix D: The soft spin Landau-Ginzburg action 


The critical action for the triangular lattice can be derived 
following the work of Blanckstein et. al.^^ Taking the unit 
vectors of the triangular lattice 

a'i = [l,0]; a '2 = [l,\/3]/2; (Dl) 

and the reciprocal lattice vectors: 

G'l = 27r[l, 1/V3]; G '2 = 27r[0, 2/V3], (D2) 

the soft modes occur at 


K± = ±[47r/3,0] . (D3) 

The soft mode expansion of the spin around these momenta is 
then given by Eq. 15. The generators of the symmetries (of the 
Hamiltonian) on the triangular lattice are: (1) unit translation 
along ai, (2) unit translation along a 2 , (3) global Z 2 . 

The two amplitudes transform as: 
















13 


% : {x,y} ^{x + l,y}-. y^v) ^ M^(r + ai) : ^ {V’+,V’-} ^ (D4) 

7y : {x,y} ^ {x,y + 1} : //"(r) ^ y^{v + SL 2 ) ■ ^ {V’+,V’-} ^ (D5) 

Z 2 : {x, 2 /} ^ {x,y} : n%r) -y\r) : ^ {V’+,V’-} (D6) 

Ce : {a:;,!/} -)> {-y,x + y} : //^(r) ^ {V’+,V’-} {V’-,'*/’+} (D7) 

J : {a;, y) {-x, -y} : /u^(r) ^ {V’+, V’-} {V’-, ^+} (D8) 


The resultant Euclidean Landau-Ginzburg action is given by 
Eq. 16. 

Eor our case we have two copies of the triangular lattice, 
which comprise the honeycomb lattice, and the two triangu¬ 
lar lattices transform into each other under inversion about the 


midpoint of the bond joining the sublattices of the medial hon¬ 
eycomb lattice. Also the inversion about the bond (J) is no 
longer there and the Cq symmetry is replaced by a C 3 sym¬ 
metry and finally there is refiection about the vertical bond 
(IR). Under these new symmetries the transformation of the 
soft modes are: 


C 3 : {x,y} -)> {-x-y,x} : : => } -)> {'ip+\'ip‘^y}; 

3b : {x,y} {-x,-y} : M^(l,r) -> y^2,-r) : ^ -)> 

: {x,y} {-x -y,y}: M^(l,r) ^ {ip^+\y^y} -5> (^9) 


This gives the action given by Eq. 19. Appendix E: Classical action on stacked honeycomb lattice 

The Hamiltonian in Eq. 14 has the following form: 

H = (El) 

{{ik)) i (ij) 

To obtain the Euclidean action we resolve the partition 
function Z = in the time direction 

by inserting compete basis of states at each time slice 

Z = ^ ^ ^ ^ Mo,l • • • Mr,r+1 • • • '^N,0 

{f^o} 

(E2) 


where 


To simplify this we first rewrite 


M, 


r,r+l ^ 


-Ara 


{{ik)) 


ArAE(,,)^r/^-eArrE,Mr 


IK+i}) 




nu 

(ij) V 








IK+i}) 


oc e 






,-2ATrEi 




■^IK+i}>- 
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1 — 1 — 


where tanh(ArA) = e Ar = P/N. Inserting a basis of and using the identity 2 " 2 " , this becomes 

“+^11 




-2Arry^. 2^2+ 2 


Gathering terms as 

{’T"ij,'r} {/^rl 


the sum over {/i^} gives 




l-^f,T I l-^f,T + l 


+Xlj 




Introducing link variables r]ij^r = ±1 on space-like links through riij^r = —this becomes 


(E3) 


oc J2 11 1^ ' 

{Vij,T} * 


^-2Arr 




r+1 


Now writing tanh(Ar) = e we get 


Mr^r+i OC ^ ]^ I cosh(^) + sinh(^) 


{'nij,T} 


Vij,T j 


'i,r+l 


OC e 


-AraJ2((i 


{{ik)) 


^ir^kr ^ ^ ^ j, r+-^ ( O j ’Vij, r ) ^ Mi, t--|-1 


{^zjjT } 

Finally gathering all terms we write the partition function as 


Zoc ^ ^ e 


-JE 


where E is the classical Hamiltonian which takes the form 


(E4) 


(E5) 


E= Y1 ^^trEkr - f E ( n ) Ei,r^lr+l “ ^ E 


{{ik)),r 


J 




where J = Ara. 


(E6) 
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